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Abstract 

This paper presents an algorithm for estimating the angular rate vector of a satellite which is based 
on the time derivatives of vector measurements expressed in a reference and body coordinate. The 
computed derivatives are fed into a special Kalman filter which yields an estimate of the spacecraft 
angular velocity. The filter, named Extended Interlaced Kalman Filter (EIKF), is an extension of the 
Interlaced Kalman Filter (IKF) presented in the literature. Like the IKF, the EIKF is a suboptimal 

Kalman filter which, although being linear, estimates the state of a nonlinear dynamic system. It 

consists of two or three parallel Kalman filters whose individual estimates are fed to one another and 
are considered as known inputs by the other parallel filter(s). The nonlinear dynamics stem from the 
nonlinear differential equation that describes the rotation of a three dimensional body. Initial 

results using simulated data, and real RXTE data indicate that the algorithm is efficient and robust. 

I. INTRODUCTION 

Small inexpensive satellites which do not carry gyroscopes on board still need to know their angular 
rate vector for attitude determination and for control loop dam ping. The same necessity exists also in 
gyro equipped satellites when performing high rate maneuvers whose angular rate is out of range of the 
on board gyros. While the attitude determination task requires high precision angular rate 

measurements, low precision angular rate measurements are adequate for control loop damping. Satellites 
usually utilize vector measurements for attitude determination. Such measurements are, for example, of 
the direction of the nadir, of the sun, of the magnetic field vector, etc. The vector measurements can 
be differentiated in time in order to obtain valuable information. This approach was used by Natanson 1 
for estimating attitude from magnetometer measurements and by Challa, Natanson, Deutschmann and Galal 2 
to obtain attitude as well as rate. 

Angular rate can be extracted from vector measurements in the following way. Let b represent a vector 
measured by an attitude sensor such as Sun Sensor, Horizon Sensor, etc. For the time being let us 
assume that b is the earth magnetic field vector. From the laws of dynamics it is known that 

i b 

b = b + © x b (1) 

i b 

where b is the time derivative of b as seen by an observer in inertial coordinates (i), b is the time 
derivative of b as seen by an observer in body coordinates (b), and © is the angular rate vector of 
coordinate system b with respect to coordinate system i. (Note that the choice of the inertial 
coordinate system as the reference coordinates is arbitrary). We can write (1) as follows 
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b i 

[bx]co = b - b (2.a) 

where [bx] is the cross product matrix of the pleasured vector b. Note that b is computable since b is 

usually known from Almanac or a model, and b is computable from the measurements. Consequently, all 

elements of (2.a) other than co are known. Let us. resolve (2.a) in the body coordinates and let us also 

denote the transformation matrix from i to b by D* then (2.a) can be written as 

0 

[bx]co = b - D* (2.b) 

0 

1 J 

where the dot denotes a simple time derivative. Note that b is resolved in the i coordinates and D fe has 

to be known. We realize that co cannot be determined from (2.b) since [bx] is not invertible. If we add 
though one more vector measurement, c, from an additional sensor, then to can be determined as shown 
next. Similarly to (2.b), we can write for c . 

[cx]co = c - D* c (3) 

When we augment (2.b) and (3) into one equation we obtain 
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d = Geo (6) 

Next, define G , the pseudo-inverse of G, as follows 

G # = (G T G)' 1 G T (7) 

where T denotes the transpose, then co, the best estimate of co in the least squares sense, is given by 

oo = G # d (8) 

Note that this solution exists only if b and c are not co-linear. An estimate of co, better than that 
given in (8), can be obtained when the problem is treated as a stochastic problem and some kind of 
filtering is applied to the measurements. Moreover, filtering in the sense of estimation is a must when 
at each time point we have only one vector measurement. (Such case exists, for example, when we use a 
Sun Sensor and some other vector measuring sensor, and the satellite happens to be in a shadowed zone). 
In such case we use the vehicle dynamics for propagating the estimate of co. As will be shown in the 
ensuing, the dynamics model of a spacecraft (SC) is a non-linear model, therefore a linear Kalman 
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filter (KF) is not suitable, and some kind of non-linear estimator is needed for estimating co. The 
extended Kalman filter (EKF) is, then, the natural choice. However, Algrain and Saniie 4 introduced the 
Interlaced Kalman filter (IKF) which is a sub-optimal filter that is a combination of two linear Kalman 
filters that operate simultaneously and feed one another. While the IKF of Algrain and Saniie was an 
ingenious idea, they did not utilize its full power since they fed the filter with the angular rate 
vector itself as measured by gyros and not with vector derivative information. Therefore they 
practically used the IKF merely as a low pass filter and not as an estimator. This is equivalent to 
using the IKF for filtering co computed in (8). In contrast to Algrain and Saniie, we use their idea to 
estimate the angular rate vector directly from vector measurements and their time derivatives and are 
able to obtain estimates even when we have a single measurement at each time point. We also extend 
their dynamics model farther to include products of inertia. This leads to the use of two or three more 
sophisticated KFs that make use of three dynamics models. We call the extended filter: Extended 
Interlaced Kalman Filter (EIKF). Finally^ our work differs considerably from that of Natanson 1 , and 
Challa, Natanson, Deutschmann and Galal mainly because most of our investigation is dedicated to the 
filtering stage. In the next section we develop the dynamics models which give rise to the use of the 
EIKF. This leads to the development of measurement equations that correspond to the states of the 
dynamic models. This is done in Section DL In Section IV we present the stochastic models which are 
used by the EIKF. They are based on the dynamics and measurement models derived in Sections II and III 
respectively. Then in Section V we introduce several options for implementing the EIKF followed by test 
results of the EIKF which we show in Section VI. Finally, in Section VII we present our conclusions 
from this work. 


H. SPACECRAFT DYNAMICS 


In order to apply a recursive estimator to estimate the angular rate vector of a gyro-less spacecraft 
(SC), one needs to utilize the dynamics model of the SC. The angular dynamics of an SC is given in the 
following equation 


Ico + h + co x (Id) + h) = T 


(9) 


where I is the moment of inertia matrix, co is the angular velocity of the satellite with respect to 
inertial space, h is the angular momentum of the momentum, or reaction, wheels and T is the external 
torque applied to the SC. All vectors in (9) are resolved in the b system. Since I is nonsingular, we 
may write (9) as 


<i> = -v\(o x (ico + h)] + r'cr-h) 


( 10 ) 


Hie inertia matrix, I, is given by 0 
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( 11 ) 


where 1^, I , and 1^ are the moments of inertia about the body major axes x, y and z respectively, 
and 1^, 1^ and I are the product moment of inertia terms. Using these notations, (10) can be 
written as follows. Define 
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X = 


then (10) can be written as 
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(14.a) 

then (13) can be written as follows 
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(14.c) 


f = i^cr-h) 


(14.d) 

(15) 


® = F «>® + V + + t 

The latter is the desired rotational dynamics equation which expresses the time derivative of CO, the 
angular velocity vector of the SC with respect to inertial space, in terms of the known forcing 
function, f, and co itself. This equation is the central equation in the development of the filter. We 
realize that the solution of (15) hinges on our knowledge of % ^d X. As will be shown later, they will 
be estimated by their own estimator. Those estimators will each need a dynamics model for the the 
vector it is set to estimate. The derivation of the dynamics model is presented next. First we 
differentiate (12.d) to obtain the second dynamics equation 
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(17.a) 

then (16) can be written as 

X = F Z X + B x “ (18 > 

which is the desired equation. To obtain the dynamics equation for X, we differentiate (12.e). This 
yields 
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Let 
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F^= 0 (20.a) 


then (19) can written as 


and 


B x = 


2co 0 


0 2a 0 


0 2(0 


X = F^X + B^co 


(20.b) 


( 21 ) 


Equations (15), (18) and (21) are the deterministic dynamics equations which describe the behavior in 
time of (0 and the products of its components. They form the foundation of the stochastic dynamics model 
of the EIKF.Next we develop the measurement equations which will serve as the basis of the stochastic 
measurement model used by the EIKF to update its estimates. 


ffl. MEASUREMENT EQUATIONS 
m.l Raw Vector Measurements 

We start by deriving the measurement equation for the primary KF whose dynamics is given in (15) and 
which estimates ( 0 . Re-write (2.b) 


Let 


z = b - D b 
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(22.a) 


then (2.b) can be written as 


[bx](o = b - D‘ b 
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and 
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C b = [bx] 


(22.b) 


(22.c) 


The measurement vector z is a computable three dimensional vector which is data to be fed into the 

(Ob ^ 

EIKF part that estimates ©. We note that C is a 3x3 singular matrix. It is obvious that one of the 

b 

three equations of (23) is a linear combination of the other two and, thus, is superfluous. Although a 
white noise will be added to z^ at a later stage (see(37)) and , thus, will turn the three equations 

3 

in (23) into independent equations, the singularity of C. will be troublesome. Problems may arise in 

b 

3 3 3 pj* 3 « 

the KF, designed to estimate co, when computing its gain according to K = P C f CP C + R V . 
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C C + R I . The matrix CP C is singular, and since the elements of the noise covariance 
b CO b ffl J h 01 b o’ 

. 3~ 

matrix, R 


CD 
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are rather small, the inverse yields a matrix whose elements are very large. This in turn 

yields a very inaccurate gain matrix. All that keeps the inverted matrix from being strictly singular 
is the noise covariance matrix. The physical meaning of this ill-conditioned case is that the noise is 
the added information which causes the dependent deterministic equations to be independent. This 
information is, of course, meaningless and should not be considered. As a remedy to this 
ill-conditioned case, we eliminate one of the rows of (22.c). The question is then which row to 
eliminate. It is clear that the answer to this question hinges on the value of the components of z . 

Obviously, if the SC rotates fast about the body z axis and not at all about the other two, then it 
can be seen either from (2.b) or (22.a) that the third component of a is negligible and thus the 

third row of (22.c) should be eliminated. For the sake of the ensuing presentation we will assume that 
this is the case. This check, however, has to be performed before every measurement update of the 
filter. Having made the latter assumption, define 
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(23. a) 


and 


C = 
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(23 .b) 
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COb b 

Next we have to develop the measurement equation needed to estimate % and A,. We can choose one of two 
options which are based on entirely different approaches. According to one approach we obtain the 
needed measurements from a second differentiation of a measured direction. It is well known , and 
indeed very easily shown, that, using the notation 

i b 

w = b (24.a) and u = b (24.b) 

the following holds 

D 1 w = u + 2© x u + © x b + © x(© x b) (25) 

b 
i 

where w is resolved in the i coordinates and, as before, the dot symbolizes time differentiation 
performed in the b coordinates. Let 


then (25) can be written as 
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then it can be easily verified that (27.a) can be written as 

VV +S| v C8) 

Note that like 3 z before, 3 z , too is a computable vector which is data to be fed into the EIKF part 
© XAb 

that estimates % and A. Now, the argument that led to the reduction of the expressions in (22) to the 

corresponding expression in (23), holds here too. Consequently, we eliminate one row in the expressions 

jn (27). (As before, the row to be eliminated is determined by the relative size of the components of 

z . ). Assume that here too, the third raw is eliminated, then if we define 
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it can be easily verified that (28) can be written as 

*xto- M b * + V (2M) 

If we now have a second vector measurement, say c, then we get an identical set of equations when now b 
is replaced by c and the subscript b is replaced by the subscript c. Explicitly, we do the following. 
Define 
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then the measurement equations are 


z = C 0 ) 

0>c c 

= M x + N X 


From the above, the extension of the measurement equations in case that we have more than two vector 
measurements at one time point is obvious. The other option for obtaining measurement equations needed 
to estimate % and X is based on co, the EIKF-generated estimate of 0). We will postpone the introduction 
of this option until we present the EIKF. 

m.2 Pre-Processed Vector Measurements 

When we measure two different vectors at the same time point, then, as shown in (8), we have enough 
equations to obtain an estimate of co without resorting to a recursive estimator. Therefore we can, 
first, compute an estimate of co using (8), and then filter the estimate using the EIKF. As mentioned 
before, this is what was basically done in [6], only that there, co was obtained as an output of 
gyroscopes rather than an analytic solution based on vector measurements. Although this approach does 
not fully utilize the capabilities of a recursive estimator, for the sake of completeness, we show here 
how to formulate the measurement equation in order to apply the EIKF in this case too. Re-write (8) 

co = G # d (8) 

Let 

*<op= G ' d 

and let U denote the 3x3 identity matrix, then we can write (30) as 

u 0) (33) 

Hie last equation is the measurement equation which corresponds to the dynamics equation of (15).The 
measurement equation for estimating % and X can be either those presented in the preceding sub-section; 
namely (29) and/or (31.b), or they can be directly related to to computed in (8). The latter will be 
explained later when we introduce the suitable EIKF. 


IV. THE EIKF MODELS 


The dynamics and measurement equations presented in Section II and Section HI respectively, are 
nominal equations. In preparing the equations for use in a filtering routine, we add to them white 
noise vectors to express model uncertainties. These uncertainties stem from two sources, first, there 
are modeling errors because the equations are not the exact dynamics and measurement models, and 
second, in the sub-optimal filter that we will use, we will assume that x and X are constant in the 
propagation time-interval that we will use to propagate the estimate of ( 0 . This assumption is clearly 
wrong even though it enables us to obtain satisfactory results. The importance of the white noise added 
to the each dynamics equation is in its PSD matrix which we adjust by trial and error to obtain the 
best filter performance. Similarly, the white noise added to each of the measurement equations 
indicates the measurement accuracy expressed by the covariance matrix of the noise vectors. This 
covariance too, is adjusted in order to yield the best filter performance. Adding the white noise 
vector, n^, to the central dynamics equation in (15), yields the following main dynamics model 

(34) 

Similarly we add white noise vectors to the right hand side of (18) and (21) which become 
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Adding white noise to the measurement equations turns, respectively, (23.c), (29.d) and (33) into 



The extension of (37) and (38) to the case where we have more than one vector measurements is obvious. 
As mentioned before, several measurement models which are based on the estimate of to, will be 
introduced when we present the EIKF in the next section. 


V. THE EXTENDED INTERLACED KALMAN FILTER 

Given the models of the preceding sections, we have several options for designing an EIKF. Like the 
models, the EIKF itself can be divided into two basic categories. The first is one which handles raw 
vector measurements, and the second category is one which handles pre-processed vector measurements. In 
the ensuing we only present the models to be used in the interlaced linear KFs. The KF algorithm itself 
can be found, of course, in standard KF texts. 

V.l Raw Vector Measurements 


We have several options for designing an EIKF when given raw vector measurements. 

The following are some options. 

Option 1 : 

We run three parallel linear KFs. The equations of the filters are as follows. 

Filter^ 1_ 

The dynamics equation is derived from (34) and the measurement equation is given in (37) 

A A 

0) = F cd + Bjc + U,X + f + n 
0) (Hr 092 09 

Z = C CD + v 
COb b (Ob 

Note that % and X are inputs from the other two filters that run in parallel to Filter 1. 


(40.a) 

(40.b) 
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equations are based on a single vector measurement; namely, b. If we obtain another vector measurement, 
say c, at a certain time point, then we use (30) and (31) to generate measurement models similar to 
(40.b), (41.b) and (42.b) and perform consecutive measurement updates of the three interlaced filters, 
or we can augment the two vector measurements in each filter and perform in each of them one combined 
measurement update at that time point. The extension of this case to more than two simultaneous 
measurements is immediate. The three filter model of Option 1 is summarized in Table I. 


Dynamics 

Me a s u remen t 

• a A 

to = F to + B x+B X + f + n (40.a) 
to to* co 2 to 

X = F Z X + B^to + n^ (41. a) 

X = F^X + B^to **X (42.a) 

*<ob m C b“ > + "bib <40b) 

l X H- V “ “bit + (41b) 

‘XU' M b* “ V + T X Xb < 42 -»> 


Table I: Filter Model of Option 1 

A block diagram representation of the EIKF of Option 1 is depicted in Fig. 1. 



Fig.l : Block Diagram of the EIKF of Option 1. 


Option 2 : 

Here we run only two parallel interlaced linear KF. They are as follows. 

Filter^ J 

This filter is identical to Filter 1 of the preceding option. 

Filter 2 

In order to present Filter 2, we adopt the following definitions 
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then (41. a)) and (42.a) can be augmented into the single dynamics equation 
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and (38) can be written to suit (45.a) as 
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The EIKF model of Option 2 is summarized in Table n. A block diagram representation of the EIKF is 
depicted in Fig. 2. 
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Table II: Filter Model of Option 2 
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Fig. 2 : Block D i ag ram of the EIKF of Option 2. 

Option 3 : 

Recall that in Option 1, as well as in Option 2, we had to use the second time derivative of the 
measured vectors in order to generate the data for the measurement models. We can use a different 
approach though that does not require a second differentiation. We simply use 0) which is estimated in 
Filter 1 and treat it in the other parallel Alters as a "measurement” of % and X for they are 
functions of (D (see (12.d.e)). Doing so we obtain the following measurement equations 
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which we can write as 


and 



(47. a) 


(47 .b) 


where, as before, U is the identity matrix. We added the white noise vectors, v and v. since on the 

X A 

left hand side of the above equations we do not have x and X, but rather their estimates. In this 
Option The three parallel filters are as follows. 

Filter 1 


This filter is identical to Filter 1 in Option 1. 


Filter^ 2 

The dynamics equation of this filter is exactly like the one given in (41. a), but the measurement 
equation is that given in (47.a); that is, 

X = F X X + (4 8. a) 

» X «UX + v z (48.b) 

Filter^ 3 

The dynamics equation of this filter is exactly like the one given in (42.a), but the measurement 
equation is that given in (47.b); that is, 

X = F^X + Bjffl + (42.a) 

Zy= U X + v^ (47. b) 

The EIKF model of Option 3 is summarized in Table III. 


Dynamics 

Measurement 

(b = F to + B x + B X + f + n (4 0. a) 
to to* to 2 to 

X = * % X + B^w + (4 1 .a) 

X = F^X + B^co + (4 2 .a) 

z = C to + v (40.b) 

COb b tOb v ' 

Z z = U x + (47. a) 

z x = U X + (47. b) 


Table III: Fi 1 ter Model of Opt i on 3 
A block diagram representation of the EIKF of Option 3 is depicted in Fig. 3. 
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Fig.3 : Block Diagram of the EIKF of Option 3. 


V.2 Pre-Processed Vector Measurements 


As we have already seen, pre-processed vector measurements yield an estimate of to, and as mentioned 
earlier, the full advantage of a recursive estimator is not utilized when a measurement or an estimate 
of © is available; however, for the sake of completeness, we present an EIKF scheme for this case too. 
The filter model of this case is similar to the model of Option 3. The dynamics equation of the present 
Filter 1 is identical to that of Option 3 but the measurement equation is different. Now the 
measurements that are fed into Filter 1 are not vector measurement, but rather a preliminary estimate 

of ©, which we denote by © , thus following (8) and (32) we define 

P 

© = G # d (48.a) 
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and 



(48. b) 


and then, following (39), we write the measurement equation of Filter 1 as 

z = U © + v (48. c) 

©p © 

As for Filters 2 and 3, while their dynamics equations are identical to those of Option 3, their 
measurement equations can be based on either the input ©^ (which is also the input to the present 

Filter 1), or on © which is the input to Filters 2 and 3 of Option 3. The EIKF of this case is as 
follows. 


The EIKF for the Pre-Processed Vector Measurements: 


We run three linear filters in parallel. 


Filter 7 

The dynamics equation is identical to that of Option 3. The measurement equation is 

z = U © + v (48.c) 

©p © 

Filter^ 2 

The dynamics model is identical to that of Filter 2 of Option 3. As for the measurement model, define 
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The dynamics model is identical to that of Filter 3 of Option 3. As for the measurement model, define 
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(50.b) 


The model of the EIKF for the pre-processed vector measurements is summarized in Table IV 


Dy nam i cs 

Measu rement 

co = F<b+By+B Z+f + n (4 0. a) 

0) G)* CO 2 00 v ' 

Z = F Z Z + B^co + (4 1 . a) 

Z = F^Z + B^co + (4 2 . a) 

z = U co + v (48.c) 

COp co v 1 

V u * + v » (4W » 

%= u x < s »*» 


Table IV: Model of the EIKF for Pre-Processed 
Vector Measurements 

A block diagram representation of the latter EIKF is depicted in Fig. 4. As mentioned earlier, here too 

we have several options. We can, for example, use co in the dynamics equation of Filters 2 and 3 in 

P 

addition to using it as measurements in these two filters. 



Fig.4: Block Diagram of the EIKF for Pre-Proces sed Vector Measurements 
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Fig. 5: True RXTE Angular Velocity Components 



Fig. 6: Estimated RXTE Angular Velocity Components 


V. FILTER TESTING 

As a first step in the testing of the EIKF for estimating to, we applied the filter presented as Option 
1 (see Table I and Figure 1) to simulated data. After obtaining satisfactory results we applied the 
filter to real data obtained from the RXTE satellite which was launched on Dec. 30, 1995. We used the 
downlinked magnetometer data (b.) and Sun sensor data (c.) as well as the wheel momentum data. We 

applied the EIKF just before the beginning of a maneuver; namely, at 21h, 43min and 31.148sec on Jan. 
4, 1996. The true rates, estimated rates, and the estimation errors are shown in Figs. 5, 6, and 7, 
respectively. 

VL CONCLUSIONS 

In this paper we presented an algorithm for estimating the angular velocity of a rigid body like 
satellite. The algorithm is based on vector measurements and their derivatives. The algorithm is an 
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extension of an estimator named, Interlaced Kalman Filter (IKF), which was introduced in the past by 
Algrain and Saniie. The IKF enables the use of several linear filters running in parallel for 
estimating the state of a non-linear dynamic system. In this paper we developed an IKF for a more 
general dynamic model and named it Extended Interlaces Kalman Filter (EIKF). Unlike Algrain and Saniie, 
we make a full use of the estimator in that we use direction vectors, rather than measured angular 
velocity to obtain an estimate of the angular velocity. In this paper we presented several versions of 
the EIKF for angular velocity estimation. 

Simulation results indicate that the EIKF is an efficient and stable estimator of the angular velocity 
vector. 



Fig. 7: Estimation Error of the RXTE Angular Velocity Components 


REFERENCES 

1. Natanson, G., "A Deterministic Method for Estimating Attitude from Magnetometer Data Only, 
"Proceedings of the 43rd IFA Congress, Paper No. IFA-92-0036, Washington, DC, Aug. 28 - Sept. 5, 1992. 

2. Challa, M., Natanson, G., Deutschmann, J., and Galal, K., "A PC-Based Magnetometer-Only Attitude and 
Rate Determination System for Gyroless Spacecraft," Proceedings of the Flight Mechanics/Estimation 
Theory Symposium, NASA-Goddard Space Flight Center, May 16-18, 1995. pp. 83-96. 

3. Gelb, A., (ed.), Applied Optimal Estimation, MIT Press, Cambridge, MA, 1974, p. 19. 

4. Algrain, M. C., and Saniie, J., "Interlaced Kalman Filtering of 3-D Angular Motion Based on Euler’s 
Nonlinear Equations," IEEE Trans, on Aerosp. and Electr. Sys., Vol. 30, No. 1, 1994, pp. 175-185. 

5. Wertz, J. R., (Ed.), Spacecraft attitude dynamics and Control, Reidel Publishing Co., Dordrecht, 
Holland, 1999, p. 523. 

6. Agrawal, B. N., Design of Geosynchronous Spacecraft, Prentice-Hall, Englewood, NJ, 1999, p. 109. 

7. Meriam, Jl., and Kraige, L.G., Engineering Mechanics Volume 2, Dynamics, Wiley, New- York, NY, 1992, 
pp. 383-386. 


69 





